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Abstract 

Monte Carlo method within, so called, classical fields approximation is applied to one dimen- 

H ' 

sional weakly interacting repulsive Bose gas trapped in a harmonic potential. Equilibrium statis- 
tical properties of the condensate are calculated within a canonical ensemble. We also calculate 
experimentally relevant low order correlation functions of the whole gas. 
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I. INTRODUCTION 



Statistical properties of multiparticle quantum systems are at the heart of statistical 
physics. The best developed are both experiments and the theory of statistical properties 
of photons. In quickly developing physics of quantum degenerate dilute gases there are two 
aspects that make the problem of thermal equilibrium atom statistics notably different than 
that of photons. First, atoms collide and their interaction is necessary for the thermalization 
process - their approach to thermal equilibrium. Second, total number of atoms is strictly 
conserved and the experiments on Bose-Einstein condensation are performed with almost 
perfect isolation of the system from an exchange of particles with the outside world. More- 
over, the system often consists of as few as a couple of hundred atoms, thus it is far from a 
thermodynamic limit. Until now, fully understood is the statistics of an ideal Bose gas only 
. In a series of papers it was found that the three commonly used statistical ensembles 
(grand canonical, canonical and microcanonical) while giving identical prediction for the ex- 
tinction of the condensate with increasing temperature, differ in the predicted fluctuations 
of the number of condensed atoms. In particular these fluctuations calculated via grand 
canonical ensemble are absurdly large. It is worth noticing that this last observation was 

Q 

already made by E. Schrodinger |2j . The statistics of a weakly interacting Bose gas is still a 
challenge. Nearly all papers deal with an academic problem of a system of N atoms confined 
in a rectangular box with periodic boundary conditions. There are two important simplify- 
ing aspects of such a confinement. First, even in the presence of (repulsive) interaction, the 
condensate wave function is still the zero momentum component of the atomic field just as 
for the ideal gas. Second is a simplicity of the Bogoliubov quasiparticle excitations spectrum 
and an analytic relation between the annihilation and creation operators of quasiparticles 
and the corresponding operators for atoms [3( . Thus, in the Bogoliubov approximation, the 
total Hamiltonian may be written as: 

H B = E + J2<k)blb k (1) 

k 

or a sum of the energy of the condensate Eq and the sum of the independent energies of 
quasiparticles (created and annihilated by fej, and respectively) which in this approx- 
imation are treated as independent bosons with simple and well understood equilibrium 
thermodynamic properties. Thus ignoring the changing number of background condensed 



atoms one can compute the statistics of the excited thermal atoms just from the statistics 
of the independent quasiparticles 4|. This way one gets a reliable results for low tempera- 
tures only. At higher temperatures the Bogoliubov spectrum gets modified and it takes a 
Bogoliubov-Popov form [5j. In this approximation the spectrum depends on the number of 
condensed atoms rather than on their total number. Moreover, also E depends on iVo- So 
)oth parameters depend on the random variable. This presents a serious technical difficulty 
6|. It has been overcome in a self-consistent way in [7|. But none of these papers takes fully 
into account the higher order terms in the interaction Hamiltonian. Physically these terms 
describe the interaction between quasiparticles that leads to their instability. In a recent 
paper [8| we have shown how to calculate the statistical properties of the weakly interacting 
Bose gas retaining full value of the interaction energy. To this end we have proposed to use 
the so called classical fields approximation In this approximation all long wave-length 
degrees of freedom of the atomic field are stripped-off their operator character and are de- 
scribed by complex amplitudes. The question if practically all atoms may be accounted for 
within the classical fields has been answered in affirmative. In a recent paper 10| we have 
shown that with a proper choice of the short wave-length cut-off all statistical properties of 

n 

an ideal Bose gas may be reproduced using the classical fields. In [8| we have applied the 
classical fields to describe statistical properties of a weakly interacting bosons again confined 
in a box with the periodic boundary conditions. In this paper we extend the Monte Carlo 
method to an experimentally relevant one dimensional weakly interacting Bose gas trapped 
in a harmonic potential. The results are obtained with the help of the classical fields, but 
are accounting fully for the interaction energy. In this case there are both theoretical 11] 
and experimental results 12| for the limited phase coherence of a very cold atomic sample, 
'mown as a quasicondensate and also recent measurements of the local density fluctuations 
1^ | . In Section [TT] we formulate the problem and describe the proposed scheme based on the 



classic Metropolis algorithm. In Section II III we present statistical properties of the conden- 
sate. In Section IIVI we present the results for low order correlation functions of the whole 
gas. 
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II. THE METHOD 



We study a one dimensional, repulsive, weakly interacting Bose gas confined in a harmonic 



trap. Excellent realizations of such a system are available 
one dimensional Bose gas has a form: 



1311 Thus, our Hamiltonian of the 



H = J &(x) + ^mwV^j i>(x)dx + 

+ | j &(x)&(x)V(x)V(x), (2) 

The Hamiltonian is a sum of the single particle oscillator energy with mass m and angu- 
lar frequency u and a conventional contact interaction with the coupling constant g. A 
convenient base is provided by the eigenstates of the harmonic oscillators (j) n (x): 

^(x) = ^<p n {x)a n (3) 

n 

where a n annihilates an atom in the n-th harmonic oscillator state. The classical field 
approximation consists in replacing the creation and annihilation operators by classical 
complex amplitudes: 

On, a\ H> a n , a* n (4) 

In [h]] we have shown, that probabilistic properties of the condensate for one dimensional 
ideal Bose gas in a harmonic trap are perfectly reproduced by the classical fields approxi- 
mation provided the number of degrees of freedom is kept finite with the last retained state 
chosen as: 

where T is the absolute temperature and is the Boltzmann constant. This suggests 
the classical fields approach may be used also for the weakly interacting gas [8j. In this 
approximation the energy of a given configuration of the field is given as: 

K 

E({a i }) = huY,n\a n \ 2 + E int ({a i }) (6) 

n=0 

where Ei nt ({ctj}) is the quartic polynomial in the amplitudes a. Hence, the probability 
distribution of a given configuration of the classical field according to the canonical ensemble 
is: 



where Z(N, T) is the classical partition function for iV atoms at temperature T. We need to 
generate this probability distribution with a set of the amplitudes subject to the constraint 
given by the fixed number of particles: 

K 

j2k\ 2 = n (8) 

n=0 

The best known Monte Carlo realization of this distribution is given by the Metropolis 



algorithm [l4j . Before we turn to the results there are two important remarks: For in- 
teracting atoms the condensate wave function is no longer the empty harmonic potential 
ground state, thus its occupation is not |ao| 2 . Instead, following [151 ]. the identification of 



the condensate requires a diagonalization of the single particle density matrix, which in our 
harmonic oscillator representation is: 



ia j ) = TT,X n ft(n)P j {n) (9) 



Pi,j = \ a 

where the mean value is taken with the probability distribution (JTj) and the eigenvector 
corresponding to the dominant eigenvalue is the condensate. The other remark concerns the 
cut-off condition (J5J). For the interacting gas we need a higher cut-off. We know that the 
repulsive gas at zero temperature is broader then the ground state of the harmonic potential 
thus its wave function needs higher energy terms. Remembering that in the Bogolubov 
approximation the quasiparticle excitations have energies counted not from zero but from 
the chemical potential \i we postulate the modification of the cut-off condition in the form: 

Khu = fi + k B T (10) 

which will be independently verified in the next Section. 

III. STATISTICAL PROPERTIES OF ID CONDENSATE 

Equipped with the numerical scheme described in some detail in the preceding Section we 
turn now to the discussion of results. We begin with statistical properties of the condensate. 



The only earlier results were obtained with the help of the time dependent methods 
In this paper we present the first results obtained with the methods of equilibrium statistical 
mechanics. In our classical fields approximation, however, effects of quantum fluctuations, 
such as quantum depletion are missing. Throughout this paper we use the oscillator units of 
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FIG. 1. (Color online)Zero temperature spatial distribution of the ground state for various values 
of the parameter g. Comparison between the classical field approximation (data from Metropolis 
algorithm) and Gross-Pitaevski equation (from imaginary time evolution). The red solid line 
corresponds to exact result for non- interacting gas. 



position, energy and temperature: 



ftw, respectively. All calculations are performed 



' k 



for 500 atoms. Our main parameter is just the coupling constant g - all simulations are 
done for g = 0.02 (what corresponds to the 87 Rb in our units) and g = 1. As a reference 
point we always have the ideal g clS CclSC. For the ideal gas we have an exact probability 
distribution of the number of uncondensed atoms N ex (thus with remaining Nq = N — N ex 
in the condensate): 

1 N 

p(N ex , T ) = - w - n (i-e), (ii) 

? l=N ex +l 

where £ = exp (—hu/kBT). We also have its best classical fields approximation: 

Pa (N m T) = (^ T -yj (12) 

with the cut-off parameter K chosen according to (|SJ) and the Monte Carlo representation 
of the classical distribution. The detailed derivation of (j!2p is given in the Appendix. 

From these, coinciding reference points we are then departing to the largely unchartered 
territory of the interacting gas. In Fig. [T]we plot the zero temperature spatial distribution 
of the condensate for several values of the coupling. Note the standard broadening of the 
condensate wave function. At this point our method merely seeks the minimum of the 
classical energy functional (j7]) satisfying the constraint (jSJ). In fact this is nearly the same as 
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FIG. 2. (Color online) Relative average number of atoms in the ground state as a function of 
temperature T. The red solid line corresponds to the exact result for non-interacting gas, points - 
results of classical field approximation. 



the ground state of the Gross-Pitaevskii equation conveniently computed by the imaginary 
time propagation. Next we look at the mean number of condensed atoms as a function of 
temperature, Fig[2J In the ID case this number decreases gradually even in iV — > oo limit. 
As it should be, the exact ideal gas result (solid, red online) is reproduced very well by the 
corresponding classical field Monte Carlo results (green crosses). Note that for a stronger 
coupling the depletion of the condensate with growing temperature becomes very rapid. 
This plot provides a direct test of our proposed modification of the cut-off condition for the 
interacting gas. The shift of the cut-off energy by the chemical potential is the smallest 
increase that guaranties ^ > 1 and in the same time gives correct zero temperature 

T— >0 

condensate wave function. Next let us look at fluctuations. It is the fluctuations that were 
ensemble dependent for the ideal gas. In Fig. [3] we plot the variance of the condensate 
occupation probability distribution for the interacting gas again compared to the ideal gas 
canonical result (solid, red online). We see the shift of the curve towards lower temperatures. 
As for all finite systems, it is not easy to define a characteristic cross-over temperature based 
on the depletion plots as in Fig. [2J It is natural to define such a cross-over characteristic 
temperature as the one corresponding to the maximum variance [3]. Of course our Monte 
Carlo method gives also access to the full probability distribution of the occupation of the 
condensate state. For T = 20 this is shown in Fig. [H 
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FIG. 3. (Color online) Relative fluctuation of number of atoms in the ground state. The red 
solid line corresponds to exact result for non-interacting gas, points - results of classical field 
approximation . 
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FIG. 4. (Color online) Probability distribution of the excited states' occupation. Solid red line - 
exact results for non-interacting gas, symbols - results of the classical field approximation. 

IV. LOW ORDER CORRELATION FUNCTIONS OF ID BOSE GAS 

In spite of all the theoretical effort, until now there are no experimental measurements 
of the temperature dependent statistical properties of the condensate. The relative fluctu- 
ations are significant only for small samples of several hundred atoms. Sufficiently precise 
determination of the condensate fraction with nearly perfect control of the total number of 
trapped atoms remains a challenge. However, measured indirectly were temperature depen- 
dent coherence properties of nearly ID Bose gas [l9( and also local density fluctuations of 
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such a gas [13j. Those are of course, also accessible to us. There are two remarks in order: 
First, our cut-off was optimized to reproduce statistics of the condensate. This way it is 
not optimized for the thermal atoms. In fact, as was shown in [h]], it overestimates the 
population of thermal modes. It must be so since in the classical field approximation all 
atoms are distributed over a finite number of low-lying states instead of the infinite number 
of them. Thus, particularly at high temperatures, when most atoms are thermal, we expect 
a deteriorating accuracy of the method. Second is the problem of ordering. Of course within 
the classical fields approximation fields do commute. Thus, there is no difference between 
the density-density correlation function and the fourth order normally ordered correlation 
function: 

But at very high temperatures the last term, missing for classical fields, becomes dominant. 
It gives rise to what experimenters call the shot noise. Classical fields cannot possibly 
reproduce this contribution unless it would be introduced by hand. The most interesting 
property of the lowest order correlation function for the ID Bosc gas has been first noted 



in Uj. Unlike in 3D, the coherence length in ID gas may be shorter than the length of the 



condensate. This phenomenon is called a quasicondensate. In the experiment [12| it has 
been observed indirectly: Upon expansion, quasicondensate's phase fluctuations turn into 
density fluctuations that can be directly observed in standard absorptive imaging. In Fig. 
Owe show the comparison of the condensate wave- function and the lowest order correlation 
function 

gi(—x,x) = 5 

(l^)| 2 > 

in the quasicondensate regime. More quantitative analysis is shown in Fig. [61 Here we vary 
the temperature confronting the length of the condensate with the correlation length for our 
two standard coupling constants. Note the crossing of the two curves defining the onset of 
the quasicondensate. We add that the ID ideal gas has no quasicondensation. 



Finally, we present results 



in Hanover 19] and in Orsay 



or the local density fluctuations. They were measured recently 



131 ] . In Fig. [7] we present temperature dependence of the local 



density fluctuations at the center of the trap defined as 

(5 2 n) _ (|*(0)| 4 )-(|*(0)| 2 ) 
(n) 2 (|^(0)| 2 ) 2 
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FIG. 5. (Color online)The first order correlation function and the density of the ground state. For 
easier comparison the quasicondensate density profile has been rescaled by the values at the center 
of the trap. 
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FIG. 6. (Color online)The correlation length and the size of the ground state wave function for 
two different values of interaction. Both lengths defined as full width at half maximum. 



Unlike in the statistical properties of the condensate itself, even for the ideal gas there is 
certain difference between the exact (solid, red online) and the Monte Carlo classical fields. 
As mentioned above it is a result of the distorted occupation of the thermal modes in the 
classical fields approach. The interaction reduces the density fluctuations as predicted in [20]. 
In Ref. jl3| the density fluctuations were measured in situ across a nearly one-dimensional 
sample versus the local density. The comparison of the ideal gas and the interacting gas of 
the same dependence is presented in Fig. [HJ It has the main qualitative properties just like in 
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FIG. 7. (Color online) Normalized fluctuations of the atoms' number at the center of the trap as a 
function of the temperature. The red solid line corresponds to exact result for non- interacting gas, 
points - results of classical field approximation. Note imperfect agreement of the classical fields 
for an ideal gas. 
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FIG. 8. (Color online)Local density fluctuations versus local density at the temperature T = 20. 
As the strength of interaction increases, the bosonic gas enters the quasicondensate regime. 



the experiment again showing the reduction of the fluctuations due to the interactions. The 
quantitative comparison is hard since experimental parameters of the number of particles 



or the temperature are known with large error bars. The authors of 13( explored also a 
very high temperature regime with the density fluctuations dominated by the shot noise. 
This is a regime not accessible to classical fields for two related reasons: in this regime 
all states have a very small, fractional occupation and also classical fields are missing the 
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quadratic term resulting from the commutator of the Bose fields, term essential in the high 
temperature regime. 

V. CONCLUSIONS 

In this paper we have demonstrated a power of the novel computational tool to tackle 
the equilibrium thermodynamics of weakly interacting gas consisting of Bose particles. The 
method, based on classical fields approximation to the full quantum theory exploits classical 
Monte Carlo technique to explore the phase space of the underlying finite dimensional clas- 
sical system, for which the set of canonical variables replaces the annihilation and creation 
operators. 

We have applied the method to the ID repulsive gas trapped in a harmonic potential 
calculating not only the statistical properties of the condensate but also the low-order cor- 
relation functions of the whole gas that may be directly confronted with the measurements. 
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Appendix: Partition function for the ideal Bose gas in a harmonic trap — classical 
fields 

The quantum partition function in the canonical ensemble may be written as follows: 

oo oo oo 

Z(N,fi)= J2 E---E--- e "^ E ^° fe ^E fe n fc ^ (A.l) 

no=0 ni=0 7ifc=0 

where /3 = l/fcgT and u is a frequency of the harmonic trap. The Kronecker delta in Eq. 
IA. II enforces in the conservation of the total number of particles. In the classical version the 
sums have to be replaced with integrals Y^=o ^ / d 2 afc and the Kronecker's delta turns 
into the Dirac's delta <^;£L n fc ,7V ^ {Yl^o \ a n\ 2 — N). Thus classical counterpart of the 
above formula reads: 
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2_ ,j2. / « 



Z(N,P) = J^...J ^ e -^^ 2 8 \^\a n \ 2 - N 



[]A e ^j|^ W '_jVl (A.2) 



7T 

j=0 \n=0 



where we have also taken into account the cut-off if, defined in (|SJ). When we use the 
integral form of Dirac's delta 6(x) = j- d.77 exp (—ii]x), the classical partition function 
takes the form 

Z(N,(3) = — / d77 e~ iNr > TT / ^ e -^^)l«il a . ( A .3) 

In the above we can easily compute all integrals over amplitudes ot^ 

1 f°° K 1 

Z(N,P) = - d ,e-««Jl—- _, (A.4) 



and then the appropriate contour integral over 77 which has if + 1 poles: 
where f = e - ^". 

Note that the expression Ylk^j J^k 1S J us ^ ec L ua l t° ("^) ; an d the whole sums of 
products can be rewritten in the simpler form: 

1 K 



kjrS^O (A ' 6) 



if!(/?M , 

Using the Newton formula we can further simplify the partition function, getting finally: 

Z(N,P) = — f-szM • (A.7) 



if! V. 

In the same manner we can compute the partition function for excited atoms: 



1 /l — £ Nex ^ K ~ 1 



(K-iy. v ^fej 

where N ex is the number of atoms in all excited states. Thus probability of finding exactly 
N ex excited atoms in the sample of N atoms at the temperature T = if is equal to 

x Z ex (N ex ,T) 1 f\-^ ex \K-i 
P d (N ex ,T)-- ' ' 
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